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Abstract 

Dipole and quadrupole solitons in a two-dimensional optically induced defocusing 
photonic lattice are theoretically predicted and experimentally observed. It is shown 
that in-phase nearest-neighbor and out-of-phase next-nearest-neighbor dipoles exist 
and can be stable in the intermediate intensity regime. There are also different 
types of dipoles that are always unstable. In-phase nearest-neighbor quadrupoles are 
also numerically obtained, and may also be linearly stable. Out-of-phase, nearest- 
neighbor quadrupoles are found to be typically unstable. These numerical results are 
found to be aligned with the main predictions obtained analytically in the discrete 
nonlinear Schrodinger model. Finally, experimental results are presented for both 
dipole and quadrupole structures, indicating that self-trapping of such structures 
in the defocusing lattice can be realized for the length of the nonlinear crystal (10 
mm). 
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1 Introduction 



Self-trapping of light in photonic lattices optically induced in nonlinear pho- 
torefractive crystals, such as strontium barium niobate (SBN) has attracted a 
considerable amount of attention, ever since its theoretical inception [1] and 
experimental realization [2,3,4]. This is perhaps primarily because it consti- 
tutes a setting where it is very natural to consider the competing effects of 
nonlinearity with diffraction and to study the effects of periodic "potentials" 
on solitary waves; here the role of the "potential" is played by the ordinary 
polarization of light forming a waveguide array in which the extra-ordinarily 
polarized probe beam evolves. 

This setting has provided a fertile ground for the experimental realization and 
detailed examination of many interesting concepts of nonlinear wave physics, 
including, for instance, the formation of discrete dipole [5], necklace [6] solitons 
and even stripe patterns [7], rotary solitons [8], discrete vortices [9] or the 
realization of photonic quasicrystals [10] and Anderson localization [11]. It is 
clear from these findings that this setting can serve not only as a host for 
the examination of localized structures that may be usable as carriers and 
conduits for data transmission and processing in all-optical communication 
schemes; it is also relevant as an experimentally tunable playground where 
numerous fundamental issues of solitons and nonlinear waves can be explored. 

It is interesting to mention in passing that in parallel (and nearly concur- 
rently) to this direction of photorefractive crystal lattices, a number of other 
venues were developed in optical and atomic physics, where the interplay of 
nonlinearity with periodicity is important for the observed dynamics. Such 
contexts involve on the optical end, the numerous developments on the ex- 
perimental and theoretical investigation of optical waveguide arrays; see e.g. 
[12,13] for relevant reviews. On the atomic physics end, the confinement of 
dilute alkali vapors in optical lattice potentials [14] has similarly offered the 
opportunity to examine many fundamental phenomena involving spatial peri- 
odicity, including the manifestation of modulational instabilities, Bloch oscil- 
lations, Landau-Zener tunneling and gap solitons among others; see [15] for a 
recent review. 

In the present study, motivated by SBN crystals and photorefractive lattices, 
we focus on two-dimensional periodic, nonlinear media. In such settings, the 
vast majority of studies has centered around focusing nonlinearities; the rea- 
son for this is two-fold: on the theoretical end, the topic of focusing solitary 
waves is interesting due to the collapse (and concomitant arrest by the lat- 
tice) phenomena in the case of cubic nonlinearity [16]. On the other hand, 
in the photorefractive case, collapse is absent due to the saturable nature of 
the nonlinearity; however, in the latter case, it is technically easier to work 
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with voltages that are in the regime of focusing rather than in that of the 
defocusing nonhnearity (in the latter case, sufficiently large voltage, which is 
tantamount to large nonhnearity, may cause damage to the crystal). Hence, 
the only coherent structure that appears to be explored in the defocusing 
regime experimentally appears to be that of fundamental gap solitons excited 
in the vicinity of the edge of the first Brillouin zone [2]. More complex gap 
structures appear not to have been studied systematically, perhaps partially 
due to the above reasons. In an earlier work [17], motivated by the above 
lack of results, we considered more complicated multipole (i.e., dipole and 
quadrupol^) and vortex structures in the defocusing case for a standard dis- 
crete dynamical lattice (namely, the discrete nonlinear Schrodinger equation 
with cubic nonhnearity [18]). This provided us with an analytically tractable 
and numerically accessible roadmap for the study of some of the relevant solu- 
tions. In the present work, we employ a continuum model involving a periodic 
potential and a saturable nonhnearity as associated with the SBN crystals. 
In particular, we demonstrate both experimentally and numerically multipole 
(dipoles and quadrupoles) solitons in 2D square lattices induced with a self- 
defocusing nonhnearity. We numerically analyze both the existence and the 
stability of these structures and follow their dynamics, in the cases where we 
find them to be unstable. We also qualitatively compare our findings with the 
roadmap provided by the discrete model [17]. 

Our presentation is structured as follows. In section 2, we present our theoret- 
ical setup (related to our experiments). Dipole solutions with the two excited 
sites in adjacent wells of the periodic potential (nearest-neighbor dipoles) are 
studied in section 3. Subsequently, we do the same for next-nearest-neighbor 
dipoles (excited in two diagonal sites of the square lattice) in section 4. Section 
5 addresses the case of quadrupoles, and section 6 presents our experimental 
results for the structures examined. Finally, in section 7, we summarize our 
findings, posing some interesting questions for future study. 



2 Setup 



For our theoretical considerations, we use the non-dimensionalized version of 
the photorefractive model with the saturable nonhnearity, as developed in 



It should be pointed out here that although in the standard electromagnetic 
convention, one thinks of dipoles as bearing two opposite charges, and quadrupoles 
as bearing two positive and two negative charges (in our case, phases of the excited 
sites), here we will use a different notational convention. More specifically, following 
the etymology of the words, we will refer to any two-site excitation as a dipole, and 
to any four-site excitation as a quadrupole herein. 
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detail in [19,20], in the following form: 



Eq 



u. 



(1) 



l + Iol + \U\^ 



In the above expression U is the slowly varying amplitude of the probe beam 
normalized by the dark irradiance of the crystal Id^ and 



is a square optical lattice intensity function in units of 7^. Here Iq is the 
lattice peak intensity, z is the propagation distance (in units of 2fciD^/7r^), 
(x, are transverse distances (in units of D/n)^ Eq is the applied DC field 
voltage (in units of 7T'^{kln^D'^rs3)~^)^ D is the lattice spacing, ko = 27r/Ao 
is the wavenumber of the laser in the vacuum, Aq is the wavelength, rie is 
the refractive index along the extraordinary axis, ki — kon^^ and r^s is the 
electro-optic coefficient for the extraordinary polarization. In line with our 
experiment, we choose the lattice intensity /q = 5 (in units of Id). A plot of 
the optical lattice is shown in Fig. 1, also for illustrative purposes regarding 
the location where our localized pulses will be "inserted". In addition, we 
choose other physical parameters consistently with the experiment as 

D = 25 /im, Ao = 0.5 /im, rie = 2.3, r^s = 280 pm/V. 

Thus, in this paper, one x oi y unit corresponds to 7.96 /xm, one z unit 
corresponds to 3.66 mm, and one £"0 unit corresponds to 12.76 V/mm in 
physical units. In the experiments, the applied voltage is — 550V/5mm, which 
gives Eq — 8.62 in our numerical simulations. 

It should also be noted here that in our experimental results the diffraction 
length can be approximately evaluated (for the beams widths and wavelengths 
that we typically used) as being 2.5mm. As our crystal extends over 10mm 
in the z-direction, it is clear that the patterns that we observe are over a few 
(roughly 4) diffraction lengths and hence if they are self-supported within such 
length scales, this will indicate that they are indeed self-trapped beams. On the 
other hand, as we will see in what follows (in our numerical simulations), for all 
the configurations that we will find to be unstable, the instability development 
will arise typically for dimensionless propagation distances of 10 < 2; < 100. 
Since these distances are considerably longer (in dimensional units) than the 
propagation distance in our crystal, all of the patterns presented in sections 3- 
5 below (even the most unstable ones) should, in principle, be experimentally 
observable in our setting. This is corroborated by our experimental results in 
section 6. 




(2) 
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The numerical simulations are done with a uniform spatial mesh with Ax = 
1/3 and domain size 30 x 30, i.e. 91 x 91 grid points for most configurations 
(see Fig. 1 for a schematic of the spatial configurations). For some of the 
configurations, a larger domain was required, in which case a domain of size 
60 X 60, i.e. 181 x 181 grid points, was used. Regarding the typical dynamics 
of a soliton when it is unstable, we simulate the z-dependent behavior using 
a Runge-Kutta fourth-order method using a step Az = 0.00025. 




-6 -4 -2 2 4 6 

X 



Fig. 1. (Color online) A spatial (x-y) contour plot of the ordinary polarization 
standing wave [lattice beam in Eq. (1)]. The localized pulses will be sitting at the 
minima of the lattice field, as opposed to the focusing nonlinearity lattice field, 
where they reside at the maxima. Points A, C, and D are used for naming 
the dipole configurations. A is a nearest-neighbor minimum of B and D (in our 
"diagonally-oriented" lattice), while it is a next-nearest-neighbor of C. Because of 
that, dipoles whose two lobes are at A and B (or at C and D) will be called 
nearest-neighbor dipole solitons. The ones that are sitting at A and C (or at B and 
D) will be called next-nearest-neighbor dipoles. 

Stationary solutions of Eq. (1) are sought in the form of [/(x, ^, z) = u{x^ y)e^^^^ 
where /i is the propagation constant and 2i is a real valued function satisfying 

fJ'U - (Uxx + Uyy) " . 2 ^ = 



The localized states u{x^y) of (3) were obtained using the Newton-GMRES 
fixed point solver nsoli from [21] and continuation was used as a function of 
/i, to follow the relevant branches of solutions. 

The propagation constant ji we consider in this report is in the first spec- 
tral gap. Using Hill's method for the 2D problem [22], for parameter values 
mentioned above, we find the band gap to be 4.2 ^ /x < 5.4qfJ. 

^ We note that the linear spectrum for ii quoted here and shown in the following 
images was the most accurate we computed using Hill's method based on consistency 
over Ax in the limit Ax 0. The linear spectrum for the discretized problem with 
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The power or the norm of the sohtary waves is defined as follows: 

DO oo "I 1/2 



J \U\^dxdy 



(4) 



We analyze the linear stability of soliton solutions u{x^ y) of (3) by applying an 
infinitesimal perturbation to them. Writing U (x, z) = e^^^ (^u{x^ y) + e^^u{x^ y 
the perturbation u{x^ y) will then satisfy the following linearized equation 



En 



(1 + lol) U — U^U* 



(5) 



where the superscript * denotes complex conjugation. We solve the above lin- 
ear eigenvalue problem using MATLAB's standard eigenvalue solver package. 

At this point, it is also relevant to summarize the results of the discrete NLS 
defocusing model of [17] that we will use for comparison with the findings 
below. These results are incorporated in Table 1, where the stability of all 
the possible combinations of in-phase and out-of-phase, nearest-neighbor and 
next-nearest-neighbor configurations is quantified in terms of their relevant 
eigenvalues of linearization. The configurations are dubbed unstable when they 
possess (for all parameter values) real eigenvalue pairs, whereas they are con- 
sidered marginally stable, when they do not always have such pairs. However, 
the latter configurations typically, in this setting, possess imaginary eigenval- 
ues with negative Krein signature (see e.g. [23] and references therein), which 
practically means that if these collide with other eigenvalues, as /i is varied, 
complex eigenvalue quartets will emerge out of Hamiltonian-Hopf bifurcations 
[24], destabilizing the relevant solution. Hence, such solutions are not always 
linearly unstable, but may become unstable for some parameter ranges. 



3 Nearest-neighbor Dipole Solitons 



In this section, we report dipole solitons where the two lobes of the wave are 
located in two nearest-neighbor (NN) lattice sites in the 2D square lattice 
potential shown in Fig. 1. The lobes can have the same phase or tt phase 
difference and are, accordingly, hereafter termed in-phase (IP) dipoles and 

our chosen discretization is slightly different however, and in particular, the first 
band edge occurs at a slightly larger value of fi. It is interesting to note that all 
saddle-node bifurcations occur at the accurate band edge, while those solutions 
which degenerate to linear Bloch modes actually degenerate beyond this value at 
the band edge particular to the discretization (not shown). 
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TP Dinolo 


Stable (N~ — ^) 


Unstable ( N — ^) 


OOP Dipole 


Unstable (TVr = 1) 


Stable {Nr ^1) 


IP Quadrupole 


Stable {Nr = 3) 


Unstable {Nr = 3) 


OOP Quadrupole 


Unstable {Nr = 3) 


Stable (A^- = 3) 



Table 1 

Summary of the stability results of the discrete model with cubic nonlinearity, stud- 
ied in [17], for all the in-phase (IP) and out-of-phase (OOP), nearest-neighbor (NN), 
as well as next-nearest-neighbor (NNN) configurations. Nr denotes real eigenvalue 
pairs and N[~ denotes imaginary eigenvalue pairs with negative Krein signature. 

out-of-phase (OOP) dipoles, respectively. Notice that due to the (diagonal) 
nature of our lattice, the nearest-neighbor configurations that we consider are 
''built" along the diagonal direction. 

3.1 In-phase Nearest-neighbor Dipole Solitons 

We have found IP dipoles in a large interval of propagation constants /i for 
the given voltage E^. We found that the solitons exist for /i smaller than 5.46, 
or for peak intensities larger than 0.144. We note that the intensity of the 
dipoles cannot be arbitrary low, a result similar to the observed results of the 
focusing case [5,20]. The relevant findings are summarized in Fig. 2. 

The top left panel of Fig. 2 shows the stability of the dipoles as a func- 
tion of the propagation constant /i, by illustrating the maximal growth rate 
(maximum real part over all perturbation eigenvalues) of perturbations. When 
max(Re(A)) = 0, this implies stability of the configuration, while the config- 
uration is unstable if max(Re(A)) 7^ in this Hamiltonian system. We found 
that the stability region of this type of dipoles is given by 4.2 < ji < 4.91, the 
left hand limit corresponding to the Bloch band. The top right panel depicts 
the peak intensity and the power of the dipoles. 

The middle left and right panels show the profile 1/ of a dipole at /i = 5 and 
the corresponding spectra at the complex plane, respectively. We see that the 
soliton is unstable due to an oscillatory instability. This is the typical insta- 
bility in this case, in line with the discrete cubic model results. Clearly, there 
is an imaginary eigenvalue with negative Krein signature [23] which upon col- 
lisions with the continuous spectrum results in Hamiltonian-Hopf bifurcations 
and concomitant oscillatory instabilities. 

As we increase ji further, the dipoles disappear in a saddle-node bifurcation. 
The bifurcation diagram is depicted in the top panels of Fig. 2. At the bifurca- 
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Fig. 2. (Color online) The top left panel shows the stability of the dipoles as a 
function of the propagation constant fi. It is stable when the spectra is purely imag- 
inary (i.e., when max(Re(A)) = 0). The top right panel depicts the peak intensity 
and the power of the dipoles. The thin line corresponds to the solution represented 
in the middle row, while the bold line corresponds to the waveform illustrated in 
the bottom row. These two branches of solutions collide and mutually annihilate 
in a saddle-node bifurcation. The shaded areas in both of these panels represent 
the bands of the periodic potential. The middle (resp. bottom) left and right pan- 
els show the profile u of the branch indicated by thin (resp. bold) line in the top 
row at /i = 5 and the corresponding complex spectral plane (Re(A), Im(A)) of the 
eigenvalues A = Re(A) + zlm(A). 
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Fig. 3. (Color online) The evolution of the dipoles shown in Fig. 2 perturbed by a 
random noise of maximum intensity 0.25% of the soliton peak intensity. Presented 
in the figure are the isosurfaces (left panels) of height 0.2 for the first configuration 
(top panels) and of height 0.1 for the second one (bottom panels) and their slices 
at some instances (right panels). 

tion point, dP/dji ^ oc, as /x ^ 5.46, i.e., at the boundary of the first Bloch 
band. At this point, the IP NN configuration collides with a configuration 
shown at the bottom panel of Fig. 2 (where the two nearest-neighbors -along 
the axis of the dipole- of the two populated wells become out-of-phase with 
them) and disappears in a saddle-node bifurcation. The corresponding profile 
and spectral plane for the saddle branch at /i = 5 is shown in the bottom right 
panel of the same figure, illustrating the strong instability of the latter. 

We have also simulated the dynamics of the solitary waves when they are 
unstable. In Fig. 3 we present the evolution of the unstable dipoles shown in 
Fig. 2. The dipoles are perturbed by a random noise with maximum intensity 
0.25% of the soliton peak intensity. Shown in Fig. 3 are the isosurfaces and 
the slices of the soliton along the propagation direction. 

The dynamics of the soliton shown in the middle panel of Fig. 2 is presented 
in the top panels of Fig. 3. One can see that even with that strong perturba- 
tion, at z = 100 the soliton still resembles its initial configuration. Physically, 
this corresponds to a propagation distance of approximately 366 mm. This 
means that the instability is very unlikely to be observed in the photorefrac- 
tive crystal lattice used in our experiments. For longer propagation distances, 
the oscillatory instability sets in and finally rearranges the dipole into a fun- 
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damental soliton type configuration principally centered around a single site. 



For the dipole shown in the bottom panel of Fig. 2, we present its dynamics 
as bottom panels of Fig. 3. We found that the instability is strong as pre- 
dicted above such that even after a relatively short propagation distance, the 
configuration turns into the more stable solution of the IP NN dipole branch. 
We did not present the further dynamics of the dipole, as it is similar to the 
upper panels. 



3.2 Out- of -phase Nearest-neighbor Dipole Solitons 



We have also found OOP dipoles arranged in nearest-neighboring lattice wells. 
We summarize our findings in Fig. 4 where one can see that the solitons exist 
in the whole entire region of propagation constant /i in the first Bragg gap, 
/i G (4.2,5.46). This smooth transition indicates that the OOP NN dipole 
solitons emerge out of the Bloch band waves; see e.g. [25] and [26] for a relevant 
discussion of the ID and of the 2D problem respectively, in the case of the cubic 
nonlinearity. Nonetheless, the OOP NN dipoles are typically unstable due to 
a real eigenvalue pair. Notice that this is in agreement with the prediction of 
the discrete model (as can be seen from Table 1). 

As the branch merges with the band edge, we observe an interesting feature, 
namely that the configuration resembles that of a quadrupole with a tt phase 

difference between two neighboring solitons, which we call H \ — quadrupoles 

below. This can be an indication that these structures bifurcate out of the 
Bloch band from the same bifurcation point. We elaborate this further in our 
discussion of the quadrupole structures in section 5. 



In Fig. 5 we present the instability dynamics of an OOP NN dipole soliton 
perturbed by similar random noise perturbation as in Fig. 3. This type of 
dipoles is typically more unstable than its IP counterpart, as is illustrated in 
the figure. In particular, in this example of unstable evolution even z — 10, 
the instability already manifests itself. One similarity of the instability of OOP 
NN dipoles with that of the IP NN ones is that the dipoles tend to degenerate 
to a single-site, fundamental gap soliton, which is stable in this setting. We 
note in passing that similar evolution results for dipoles and quadrupoles (but 
for short propagation distances) in the focusing case are discussed in [20] . 
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Fig. 4. (Color online) The top panels correspond to the same panels of Fig. 2 but 
for OOP NN dipole solitons. The bottom panel shows the profile u and the corre- 
sponding spectra in the complex plane of the dipoles at /x = 5.4. 





Fig. 5. (Color online) Similar to Fig. 3 but for the evolution of OOP NN dipoles. 
Shown are the isosurfaces of height 0.05 (red) and 0.015 (blue) and the contour plot 
slices at some select propagation distances. 

4 Next-nearest-neighbor Dipole Solitons 



We have also obtained dipole solutions that are not oriented along the two 
nearest-neighboring lattice wells, but rather where the two humps of the struc- 
ture are located at two adjacent next-nearest-neighboring lattice sites. These 
humps can once again have the same phase or tt phase difference between 
them. We will again use the corresponding IP and OOP designations for these 
next-nearest-neighbor (NNN) waveforms. Notice here that NNN configura- 
tions that we consider are among the "closest" next-nearest-neighbor pairs, 
i.e., with the structures being aligned horizontally. In principle, one can con- 
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sider more remote pairs of next-nearest-neighbors (e.g., along the diagonal), 
however the main qualitative stability properties discussed below would not 
change in such a case. 



4^1 In-phase Next-nearest-neighbor Dipole Solitons 
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Fig. 6. (Color online) The top panels correspond to the same diagnostics as in Fig. 
2 but for IP NNN dipole solitons. The bottom panels show the profile u and the 
corresponding spectral plane of the IP NNN dipole at /i = 5.1. 




Fig. 7. (Color online) Similar to Fig. 3 but for the evolution of the IP dipole shown 
in Fig. 6. Presented are the isosurfaces of height 0.1 (red) and 0.05 (blue) and the 
sliced contours at some selected propagation distances. 

We have obtained this type of IP, NNN dipole solitons in a wide parameter 
range. The stability, the power and peak intensity of these dipoles are shown 
in Fig. 6. These dipoles typically possess a real eigenvalue (again in line with 
the discrete cubic model prediction of Table 1). 
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If one compares these IP NNN dipole solitons with the IP NN dipole sohtons 
of Fig. 2, there are clear differences, such as the typical instability (with a real 
eigenvalue) of the former in contrast with the typical stability of the latter 
(possessing an imaginary eigenvalue of negative signature that can potentially 
become unstable upon collision with another eigenvalue). Another important 
difference is that the present IP dipoles emerge out of the Bloch band edges, 
contrary to what is the case for their IP NN counterparts where the soli- 
tons disappear through a saddle-node bifurcation. Interestingly, as the branch 
approaches the upper band-edge, the profile of an IP NNN dipole becomes 

similar to that of a H \ — quadrupole (see section 5.1 below) and of an OOP 

NN dipole, as the relevant limit is approached. 

In Fig. 7, we present the dynamical evolution of the IP NNN dipole shown in 
Fig. 6 under similar random noise perturbation as above. Here, we also see that 
the instability appears earlier than the IP NN dipoles. The two humps deform 
until they become one hump already at the propagation distance z ^ 20^ 
transforming the dipole solution into a single-hump gap soliton. 



4-2 Out- of -phase Next-nearest-neighbor Dipole Solitons 

We have also obtained OOP NNN dipole solitons. A typical profile of this 
family of solutions for fi = 5.3 is shown in Fig. 8. The power diagram of 
these solitons is presented in the top panel of Fig. 8. One important finding 
in this case is that there is a relatively wide stability region, i.e., typically 
these structures are stable (as indicated again by the comparison with the 
results of Table 1). We have found that the stability range for Eq = 8.62 
is 4.15 < /i < 5.17. This class of solutions typically suffers an oscillatory 
instability due to the presence of a single eigenvalue with negative signature 
and its collision with the continuous spectrum, as shown in the middle right 
panel of Fig. 8 (cf. once again with the discrete model prediction of Table 1). 

For this solution also, we observe that similarly to the IP NN dipoles, it 
disappears at a non-zero peak intensity, in particular for max(|iip) ^ 0.11. The 
disappearance is because of collision of this dipole with another configuration 
shown in the bottom panels of Fig. 8 in a saddle-node bifurcation. It is relevant 
to note that the point of the bifurcation is very close to the edge of the Bloch 
band, i.e., for /i ^ 5.46. 

Subsequently, we simulate the dynamics of the instability and present it in 
Fig. 9. One can see that for the OOP NNN, tunneling occurs during the prop- 
agation of the soliton along the z-direction, finally leading to localization at a 
well different from the original support of the dipolar structure. Regarding the 
counterpart solution of the other branch, the evolution along the ^-direction is 
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Fig. 8. (Color online) The top panels depict the largest real part of the critical eigen- 
value, as well as the power and the peak intensity of the OOP NNN dipole solitons. 
The middle panels show the profile u and the corresponding spectra in the complex 
plane of the dipole at = 5.3, and the bottom is the unstable saddle configuration, 
which collides with the OOP NNN profile in a saddle-node bifurcation, shown for 
/i = 5. 
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Fig. 9. (Color online) The evolutions of the OOP NNN dipole and its corresponding 
counterpart for the same propagation constant fi, i.e. = 5.3, shown in Fig. 6. 
Presented are the isosurface of height 0.1 and the slices at particular propagation 
distances where one can see clearly that the light tunnels away from its initial lattice 
wells. 

a bit different as there is no tunneling away from the original position. Yet, this 
configuration also leads to localization at a single site in its final configuration. 



5 Nearest-neighbor Quadrupole Solitons 



We now turn to the examination of quadrupolar structures with configurations 
having four lobes at four adjacent lattice wells. Such structures turn out to 
exist in a large parameter region as well. In particular, we will focus our 
discussion on the case of such nearest-neighbor quadrupoles. 

Before we proceed with the discussion of our findings, we will explain the 
notations that we use for the quadrupole configurations. Since there are four 
humps with the phase difference between two neighboring wells being either 
or TT, we have four possible configurations (up to rotational and phase symme- 
tries), i.e. H \ — , H h, + + ++, and H , where and represent 

the sign of the excitation of the main hump of i.e. whether it is positive 
(phase 0) or negative (phase tt), respectively, at the four excited sites A, B, C 
and D of Fig. 1. 
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5.1 H \ — Nearest-neighbor Quadrupole Solitons 



First, we consider the H \ — type. In Fig. 10, we present a field profile of 

this structure for /x = 5. We found that this configuration exists in the entire 
forbidden-band ji G (4.2, 5.46), i.e., OOP NN quadrupoles also bifurcate from 
the edge of the Floquet band. Nonetheless, we observe that a typical instability 
that these structures have is 3 real eigenvalue pairs (two of which are coincident 
in the bottom right panel of Fig. 10), in line with the predictions of Table 1. 




-15 -10 -5 5 10 15 -0.5 0.5 



X Re(?i) 

Fig. 10. (Color online) The maximum of the real part of the critical eigenvalue (top 

left panel) and the power and the peak intensity of the H \ — quadrupoles (top 

right panel). The bottom panels depict the profile of a H \ — NN quadrupole 

at ^ = 5 and the corresponding linearization spectral plane of eigenvalues. 

In Fig. 11, we present the evolution of the unstable H \ — quadrupole shown 

in Fig. 10. One can see that the configuration is strongly unstable resulting 
in a breakup of the structure already for z ?^ 10 to a IP NNN dipole state. 
Subsequent evolution also visits the other state associated with the quadrupole 
branch, namely the OOP NN and eventually results into a single hump gap 
soliton. 

5.2 + + ++ Nearest-neighbor Quadrupole Solitons 

+ + ++ quadrupoles have four in-phase humps at adjacent lattice sites. A 
typical example of this sort is shown in Fig. 12 at /i = 5. We have analyzed 
the stability of this configuration, finding that it is stable in a wide parametric 
range, namely for 4.09 ^ /i ^ 4.93, in line with what is predicted by the cubic 
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Fig. 11. (Color online) The evolution of the quadrupole presented in Fig. 10. Shown 
are the isosurfaces of height 0.15 (red) and 0.07 (blue) and its slices at some selected 
propagation distances. 

nonlinearity discrete model in Table 1. Similarly to the IP NN and the OOP 
NNN branch this branch disappears at a non-zero peak intensity (in this case, 
for max(|t^p) ^ 0.16). 
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Fig. 12. (Color onhne) The same as Fig. 10 but for the + + ++ quadrupoles. The 
top panels show the maximal instability growth rate of perturbations (left panel) 
and the branch optical power and maximal intensity (right panel), while the bottom 
panels show the solution profile (left) and linear stability (right panel) for ji — ^. 

The dynamics of the configuration when it is unstable has also been simu- 
lated; as a particular example, we integrate the field profile shown in Fig. 
12 in z. Even though it is unstable, we found that propagation even for 100 
units of (dimensionless) length, the profile still resembles its initial condition. 
This seems to indicate that such structures should be rather straightforward 
to observe experimentally. For sufficiently long propagation, the dipole will 
eventually be destroyed through an oscillatory instability. 
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Fig. 13. (Color online) The evolution of the quadrupole presented in Fig. 12. Shown 
are the isosurfaces of height 0.2 (red) and 0.15 (blue) and its slices at particular 
propagation distances. 

5.3 H h and H Nearest-neighbor Quadrupole Solitons 

One can also examine other types of quadrupole configurations such as H h 

or + + H — . However, one then typically finds, as may be inferred by their non- 
symmetric profile, that these configurations are either always unstable (as is 

e.g., the case with H h) or are not straightforward to obtain numerically 

in a steady state form (as is e.g., the case with H ). As a representative 

example of these asymmetric profiles, we show in the middle panels of Fig. 

14 the case of the H h profile, which is always unstable due to two real 

eigenvalue pairs (and potentially also oscillatorily unstable). 

Note that this configuration is equivalent to the superposition of two IP NN 
(or two OOP NN, etc.) dipoles. Another saddle node bifurcation happens in 
this case as the solution approaches the band edge, leading to the collision of 
this branch with a more extended saddle configuration shown in the bottom 
panels of Fig. 14. 

In this case also, observing the dynamics of the instability in Fig. 15, we have 
seen the relatively fast degeneration of the mode into a single site solitary 
wave profile. The corresponding evolution dynamics for quadrupoles in the 
focusing case (but for short propagation distances) can be found in [20]. 



6 Experimental Results 

In our experiments, we use a setup similar to that used in [27]. A partially spa- 
tially incoherent beam (of wavelength 488 nm) is produced by use of a rotat- 
ing diffuser. A negatively biased photorefractive crystal (SBN:60 6xl0x5mm^) 
provides a self-defocusing nonlinearity. An amplitude mask is also used to 
spatially modulate the otherwise uniform beam after the diffuser, in order to 
produce a 2D periodic lattice. The mask is then imaged onto the input face 
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Fig. 14. (Color online) The top four images are the same as Fig. 10 but for the 

H h quadrupoles; once again the thin lines of the top row correspond to the 

solution profile and spectral plane of the middle row, while the bold lines of the top 
row to the solution profile and spectral plane indicated in the bottom row. 

of the crystal, thus creating a pixel-like input intensity pattern with a spatial 
period of about 25 /xm. The ensuing lattice beam (represented by loi in our 
theoretical model) is diagonally oriented and ordinarily polarized, thus the 
induced waveguide arrays remain invariant during propagation. An extraordi- 
narily polarized beam splitting from the same laser output is used as the probe 
beam (this is the complex field u that we monitor in our analysis above) . The 
probe beam (intensity about 6 times weaker than that of the lattice beam) is 
sent into a Mach-Zehnder interferometer to create a dipole-like or quadrupole- 
like input pattern whose phase is controlled with the piezo-transducer (PZT) 
mirrors. The input/output intensity patterns and the k-space power spectra 
of the lattice and soliton- forming beams are monitored with CCD cameras. 



By employing a defocusing nonlinearity, the two-dimensional pixel-like inten- 
sity pattern induces a backbone waveguide lattice whose intensity minima cor- 
respond to the waveguide sites [2,28]. Launching two narrow Gaussian beams 
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Fig. 15. (Color online) The evolution of the quadrupoles presented in Fig. 14. Top 
panels show the isosurfaces of height 0.2 (red) and 0.15 (blue) and the contour at 
select propagation distances of the configuration in the middle panels of Fig. 14. 
Bottom panels are the corresponding figures for the counterpart soliton. The iso- 
surface shown is of height 0.1. Both configurations eventually give rise to single-site 
localization. 

into two nearest-neighbor or next-nearest-neighbor waveguide sites either with 
in-phase or with out-of-phase relations, the probe beams evolve into dipole- 
like gap solitons through the 10 mm SEN crystal under a proper strength of 
nonlinearity. Typical results are shown in Fig. 16, where the output intensity 
patterns, upon propagation through the crystal (first column), show two main 
bright spots. This indicates that the energy of the probe beams is mostly 
localized in the waveguides that were initially excited. The interferograms be- 
tween these output patterns and a tilted plane wave show that these two main 
spots remain in-phase or out-of-phase, maintaining their initial phase relation, 
although the secondary intensity peaks adjacent to the main ones are always 
out-of-phase with the primary spots. The spatial spectra (Fourier transform) 
of the patterns (third column) are in good agreement with those obtained from 
the Fourier transform of the corresponding theoretical solutions (last column) . 
It is important to note here that although some of these configurations such 
as the OOP NN and the IP NNN dipoles (second and third rows in Fig. 16) 
have been found to always be unstable, the limited propagation distances in- 
side the crystal (10mm which corresponds to z ^ 2.73 in our dimensionless 
units) are too short for the instability to develop to an observable degree in 
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the experiment, while it is clearly shown in the theoretical analysis of soliton 
solutions in the previous sections. 

We have also performed experiments to excite quadrupole-like gap solitons, 
with four narrow Gaussian beams launched into four adjacent waveguide sites 
(i.e., four adjacent intensity minima) for both in-phase and out-of-phase condi- 
tions. Under a proper level of defocusing nonlinearity, self-trapping is observed 
with four principal intensity peaks localized in the waveguides initially excited 
[Fig. 17(a)] for both in-phase (top) and out-of-phase (bottom) conditions. Fur- 
thermore, the measured spatial spectra [of Fig. 17(b)] for these quadrupole-like 
gap solitons are also in good agreement with the corresponding ones obtained 
from the theoretical solutions [as shown in Fig. 17(c)]. 
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Fig. 16. (Color online) Observations of dipole like gap soliton. From top to bottom: 
IP NN, OOP NN, IP NNN, and OOP NNN. Panel (a) shows the output intensity 
patterns, (b) shows the interferagrams of these patterns with a tilted plane wave, 
(c) shows the spatial (Fourier) spectra and, (d) shows those spectra from theoretical 
calculations of the corresponding solutions. 



7 Conclusions 

In this communication, we examined in detail theoretically and numerically 
the existence, stability and dynamics of multipole lattice solitons excited 
with a saturable defocusing photorefractive nonlinearity. In this connection, 
we have obtained a wide array of relevant structures, including dipoles and 
quadrupoles, examining both the different possible configurations (in-phase, 
as well as out phase profiles), as well as cases where the excited sites are at 
nearest or at next-nearest-neighboring lattices sites. 
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Fig. 17. (Color online) Observations of in-phase (top) and out-of-phase (bottom) 
quadrupole-like gap solitons. Panel (a) shows the output intensity patterns, (b) 
shows the spatial spectra, and (c) shows those spectra from theoretical calculations 
of the corresponding profiles. 

We have found configurations such as the in-phase, nearest-neighbor dipole, 
the out-of-phase, next-nearest-neighbor dipole or the in-phase nearest-neighbor 
quadrupole which are typically stable (although they may incur oscillatory in- 
stabilities) . We have also identified those solitons including e.g., the in-phase, 
next-nearest-neighbor dipole, the out-of-phase, nearest-neighbor dipole or the 
out-of-phase nearest-neighbor quadrupole which are typically unstable due to 
exponential instabilities and real eigenvalues. The stability results have been 
found to be in direct agreement with the results of the analysis of the proto- 
typical discrete model of this type, namely the discrete nonlinear Schrodinger 
equation. Furthermore, we have also identified an interesting set of bifurca- 
tions that are associated with the parametric continuation and termination of 
some of the above branches. The dynamical instabilities encountered in the 
present work have been monitored through direct integration of the relevant 
dynamical equation. In most cases, the result of evolution has been the degen- 
eration of the structure to a robust single site solitary wave, although more 
complicated oscillatory evolutions are also possible. 

Our experimental results agree with the theoretical predictions. Different dipole 
and quadrupole configurations were realized in the experiment in a photore- 
fractive crystal with limited propagation length, even the most unstable ones. 
This is simply because the growth rate of the dynamical instabilities is such 
that the instabilities will manifest themselves typically for propagation dis- 
tances on the order of 100- 1000mm, but they are not appreciable within the 
short distance comparable to the length of our crystal. 

Since the framework of defocusing equations has been studied far less exten- 
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sively than their focusing counterparts, it would be particularly interesting 
to extend the present considerations to other structures. Perhaps the most 
interesting example would be the study of discrete vortices, not only of single 
but also of multiple charge which would be an interesting endeavor both from 
a theoretical, as well as from an experimental point of view. Such structures 
could be directly connected to ones studied theoretically also in other fields, 
including e.g. weakly interacting Bose-Einstein condensates; relevant exam- 
ples of vortices of charges 1 and 2 are given by the tightly bound structures of 
Figs. 8 and 10 respectively of [29]. Additionally, the present considerations of 
defocusing nonlinearity may be extended to other settings, including the case 
of radial lattices of [8]. In the latter case, one would expect to observe modula- 
tionally stable, radial beam waveforms with or without vorticity [30,31]. Such 
studies are currently in progress and will be reported in future publications. 
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